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Abstract 

We present algorithms to solve relativistic hydrodynamics in 3+l~dimensional situ- 
ations without apparent symmetry to simplify the solution. In simulations of heavy-ion 
collisions, these numerical schemes have to deal with the physical vacuum and with 
equations of state with a first order phase transition between hadron matter and a 
quark-gluon plasma. We investigate their performance for the one-dimensional expan- 
sion of baryon-free nuclear matter into the vacuum, which is an analytically solvable 
test problem that incorporates both the aspect of the vacuum as well as that of a phase 
transition in the equation of state. The dependence of the lifetime of the mixed phase 
on the initial energy density is discussed. 
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1 Introduction 



Hydrodynamics represents (local) energy, momentum, and charge conservation |^. Be- 
cause of its simplicity it has found widespread application in studying the evolution of nuclear 
matter in heavy-ion collisions. The physical input is the nuclear matter equation of state 
(EoS). Therefore, hydrodynamics allows to study nuclear matter properties like the phase 
transition to the quark-gluon plasma (QGP) |0 in a simple, straightforward manner. 

One of the first applications of hydrodynamics to heavy-ion collisions was the study 
of multi-particle production by Fermi and Landau Simple hydrodynamical models 

predicted the occurrence of nuclear shock waves in the early seventies. Full 3+1- 
dimensional calculations for heavy-ion collisions at BEVALAC energies (-E^a'b = 0.1 — 2 
AGeV) were performed about twenty years ago (see the very detailed review of this topic 
in 0, and refs. therein). It was found that the compressional shock waves created in the 
collision lead to collective flow phenomena like sideward deflection of matter in the reaction 
plane ("side-splash" and "bounce-off") as well as azimuthal deflection out of the reaction 
plane ("squeeze-out"). 

One of the main successes of the fluid-dynamical picture was the confirmation of these 
collective flow effects by BEVALAC experiments 0. More detailed investigations of the 
flow, however, have revealed p] that viscosity might be important to explain these phe- 
nomena quantitatively. Also, systematic hydrodynamical studies were done only in the 
non-relativistic limit. Therefore, ideal and dissipative relativistic hydrodynamical investi- 
gations for BEVALAC energies are mandatory. A comparison with recent excellent triple 
differential data from the EOS-collaboration [^J would help to understand in how far the 
hydrodynamical picture is quantitatively applicable at these energies. 

3-|-l-dimensional ideal relativistic hydrodynamics was applied to study heavy-ion colli- 
sions at AGS {E^^^ ~ 11 AGeV) and CERN-SPS energies {E^^^, ~ 60 - 200 AGeV) [|l|. 
Since ideal hydrodynamics assumes that matter is in local equilibrium at every instant, col- 
liding fluid elements are forced by momentum conservation to instantaneously stop and by 
energy conservation to convert all their kinetic energy into thermal energy (compression and 
heating via shock waves). However, the longitudinal rapidity loss in individual nucleon- 
nucleon collisions is limited. Thus, immediate complete stopping is not achieved in reality 
and, for higher beam energies, it is no longer justified to treat the initial stage of the reaction 
in an ideal hydrodynamical picture. Ideal hydrodynamics might nevertheless be applicable 
in the expansion stage of the collision [|TT|, where the conditions of local thermodynamical 
equilibrium are more likely to be established. In order to describe the initial stage, however, 
one has to extend the hydrodynamical description to account for non-equilibrium effects. 

While dissipative hydrodynamics in principle provides the theoretical framework for the 
study of non-equilibrium effects, its viability as a causal theory is still under debate 0. 
Moreover, all known formulations which are rigorously derived from kinetic theory rely on 
the assumption that non-equilibrium effects can be treated as perturbations, or in other 

^Over the past years the original Landau model and also simplified versions (so-called "fireball" models) 
have been frequently applied to explain heavy-ion collision data. To discuss these approaches in greater 
detail is out of the scope of this paper. 
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words, that the momentum distribution of particles deviates only slightly from its local 
equilibrium form. 

This is not the case in the initial stage of heavy-ion collisions where incident energies 
can be much larger than the Fermi-momentum of the nucleons, and thus, the momentum 
distributions of target and projectile nucleons are strongly displaced. This has led to the 
development of the so-called multi-fluid-dynamical approach to heavy-ion collisions O, O . 
Here the projectile and target nuclei are treated as ideal fluids, but the coUisional interaction 
between them is calculated via kinetic theory. In the three-fluid model of Ref. |TB[ the 



energy and momentum loss in these collisions serve as source terms for a third fluid which is 
also assumed to be in local equilibrium. Thus, apart from its creation mechanism (and from 
rescattering processes while it overlaps in space-time with the first two fluids), the dynamical 
evolution of this fluid is similar to that in the Landau model. For the AGS to CERN-SPS 
energy range this model was shown to reproduce a variety of experimental observables, like 
meson rapidity distributions as well as meson and photon transverse momentum spectra [|14 



Of particular interest would be a comparison of the hydrodynamical flow predicted by this 
model with the corresponding experimental data at the AGS |15| . 



Finally, hydrodynamics was applied to the ultrarelativistic energy regime (RHIC and 
LHC energies, Ecm > 100 AGeV) as well. In this case, the nuclear stopping power is even 
for large nuclei not sufficient^ to stop at least a few nucleons completely and the nuclei 
will essentially pass through each other, leading to a baryon-free region around midrapidity 
?/ ~ 0. However, individual parton collisions copiously generate gluons which form a region 
of high energy density in that part of phase space. This picture is supported by event 
generators based on parton-parton interactions and has led to the proposal of the 

so-called "hot-glue scenario" |]I9|. After a (proper) time of order 1 fm, matter is supposed 
to establish local thermodynamical equilibrium and the subsequent dynamical evolution can 
be described by hydrodynamics. 

This picture was anticipated long ago by Bjorken [^] who suggested a one-dimensional 
hydrodynamic scaling solution (i.e., where the longitudinal velocity Vz = z/t) for the evolu- 
tion of matter in the central region. Over the years this picture was refined to include trans- 
verse motion |2l|] and the hadronization phase transition |2^. Due to its simplicity, 
it is sometimes even used to explain experimental observables for heavy-ion collisions at the 
SPS , although at these energies neither is there a baryon-free region at mid-rapidity p6| 
nor is the longitudinal motion likely to be of the scaling type. 

Recent investigations with the HIJING event generator have revealed that for RHIC 
energies the initial conditions for the hydrodynamic stage of the system's evolution might not 
be homogeneous on surfaces of constant proper time, as in Bjorken's picture [^, but subject 
to strong fiuctuations0. As was argued in hydrodynamics might still be applicable to 
study the subsequent evolution, but the assumptions of scaling for the longitudinal motion 
and cylindrical symmetry for the transverse motion are certainly no longer valid. 

All these more recent developments call for schemes to solve hydrodynamics, which are 
applicable in the general case, i.e., without obvious space-time symmetries that help to 



^See, however, the resuhs of ||lq| . 

■^ActuaUy, this possibihty was already pointed out by Bjorken 
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simplify the solution. Such schemes are necessarily of numerical nature. In a previous 
work [^], several numerical algorithms for relativistic hydrodynamics were presented and 



shown to perform well for analytically solvable test cases like the Riemann problem for 
matter with an ideal gas EoS. These investigations present a necessary first step in the 
development of a hydrodynamical scheme for simulations of heavy-ion collisions. Following 
this line we consider in this work the ability of the algorithms to cope with two problems 
specifically occurring in simulations of heavy-ion collisions. The first is the presence of 
vacuum and the second the fact that a phase transition in the EoS changes the character of 
the hydrodynamical solution. 

Vacuum in relativistic hydrodynamics poses the following problem: matter can at most 
propagate with the velocity of light into the vacuum, i.e., during a time At matter trav- 
els at most a distance At into the vacuum^. The equations of relativistic hydrodynamics 
respect this causality requirement since they are hyperbolic. On the other hand, any algo- 
rithm solving the finite difference form of a hyperbolic differential equation has to fulfill the 
Courant-Friedrichs-Lewy (CFL) criterion [^, namely that the ratio of time step width to 
spatial grid spacing is less than unity. At/ Ax = A < 1. Causal matter transport covers at 
most a distance At = AAx < Ax in one time step. However, the algorithms average trans- 
ported quantities over a cell after each time step. Therefore, even if the algorithm transports 
a certain amount of matter causally into the vacuum during a time step, it will be evenly 
distributed over a spatial distance Ax > XAx at the end of the time step. Thus, a part 
of matter is acausally propagated a distance (1 — A) Ax into the space-like light cone. We 
will henceforth call this purely numerical phenomenon prediffusion. Obviously, the larger A, 
the smaller this prediffusion. However, some Flux Corrected Transport (FCT) algorithms, 
like the SHASTA [0, require A < 1/2 and thus might lead to a considerable amount of 
prediffusion. We will consider the prediffusion of this algorithm and also of the relativistic 
Harten-Lax-van Leer-Einfeldt (HLLE) scheme [^], a Godunov-type algorithm for the 
one-dimensional expansion of matter into vacuum. 

The second problem is related to one of the main advantages of the hydrodynamical 
approach over microscopic models, namely the ability to dynamically study the phase tran- 
sition to the QGP in a relatively straightforward way: this phase transition modifies to 
first approximation (neglecting non-equilibrium phenomena such as supercooling and bub- 
ble formation [32|) only the EoS which is given input to the hydrodynamical equations of 
motion, while a microscopic model would have to implement a microscopic mechanism for 
hadronization and deconfinement ||33| . 

At present, the order of the phase transition is not clear MM. In the absence of more 
fundamental theoretical information, especially for baryon-rich matter, one commonly resorts 
to the MIT bag EoS for the QGP which is matched to a state-of-the-art hadron matter 
EoS 1^ via Gibbs conditions of phase equilibrium (see also [0). Thus, the transition is 
first order by construction. 

However, matter with a first order phase transition may exhibit thermodynamically 
anomalous behaviour in a certain range of independent thermodynamic variables. Such 



*We use natural units fi — c ^ = 1- 
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behaviour is signalled by a change of sign of the quantity 



+ 2 c- 




(1) 



where = dp/de\cr is the velocity of sound, e,p, and a are the energy density, pressure, and 
specific entropy, respectively [Q. For thermodynamically normal (TN) matter, S > 0, for 
so-called thermodynamically anomalous (TA) matter S < [^. For the expansion of TN 
matter a simple rarefaction wave is the stable hydrodynamical solution and for TA matter 
it is a rarefaction shock wave (see Section 4). On the other hand, for the compression of TN 
matter a compressional shock wave forms the stable hydrodynamic solution while for TA 
matter this is a compressional simple wave (the compression of TA matter is investigated in 
a sequel to this paper [p7| ). 

An example for matter with TA behaviour is (net) baryon-free nuclear matter with 
an EoS of a (massless) pion gas at low temperatures and a QGP (described by the MIT 
bag EoS) at high temperatures with a first order phase transition between these phases at 
the temperature T^. While the pure QGP and pionic phase are TN, the first order phase 
transition between them gives rise to a mixed phase where S changes its sign and, thus, 
where the EoS becomes T7\0. For the one-dimensional expansion of matter into the vacuum 
we investigate how well the numerical algorithms are able to accomodate to this change in 
the thermodynamic behaviour and reproduce the correct pattern of rarefaction, i.e., a simple 
rarefaction wave in the TN region and a rarefaction shock wave in the TA region of the EoS. 

Our considerations focus first on the one-dimensional expansion of semi-infinite matter, 
because in this case an analytical solution exists both for TN and TA matter, against which 
the results produced by the numerical algorithms can be compared. We then study the 
one-dimensional expansion of matter extending over a spatially finite region such as occurs 
in heavy-ion collisions and discuss the lifetime of the mixed phase as a function of the initial 
energy density. 

This paper is organized as follows. In Section 2 we review the equations of ideal relativistic 
hydrodynamics. Since we want to keep this presentation as self-contained as possible, to 
guarantee reproducibility of our results, and to facilitate and encourage the use of the two 
algorithms considered here. Section 3 gives an exact description of the SHASTA and the 
relativistic HLLE as employed here. In Section 4 we construct the analytical solutions for the 
one-dimensional expansion of both TN and TA semi-infinite matter into the vacuum. Section 
5 contains our results on the reproduction of these results by the numerical algorithms. In 
Section 6 we discuss finite expanding systems and Section 7 closes this work with a summary 
of our results. An Appendix presents modifications of the SHASTA and their performance 
for some of the test problems of Section 5. 



^In a strict sense, E only becomes zero for this EoS. For the expansion problem studied here, this 
nevertheless leads to a discontinuity in the hydrodynamic expansion, see Section 4 for details. 
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2 Ideal relativistic hydrodynamics 



The relativistic hydrodynamical equations represent (local) energy-momentum conservation 

d^.T^"' = (2) 

and (local) charge conservation 

d^N'' = . (3) 

For heavy-ion collisions, the conserved charge is e.g. the (net) baryon number, the (net) 
strangeness etc. and there is an equation of the type (^ for each conserved charge. Provided 
that matter is in local thermodynamical equilibrium, the energy-momentum tensor T^*^ and 
the charge current N'^ assume ideal fluid form [^, i.e., 

T^"' = (e + p)u>'u'' - pg^"" , (4) 
iV^ = nu^ , (5) 

where e,p, and n are energy density, pressure, and charge density in the local rest frame of 
the fluid, = 7(1, v) is the fluid 4-velocity (7 = (1 — v^)~^/^, v is the fluid 3-velocity), 
and g^'^ = diag(+, — , — , — ) is the metric tensor. The equations of ideal fluid-dynamics are 
closed by specifying an EoS for the matter under consideration in the form p = p{e, n). 

For the numerical implementation discussed in the next section it is convenient to write 
the energy density in the calculational frame as T^^ = E, the momentum density as T"* = 
M\ and the charge density as = R, respectively. The equations of motion (|^, ^ then 
take the form 

dtE + V-{Ev) = -V-(pv), (6) 
dtM + V ■ (Mv) = -Vp, (7) 
dtR + V-{Rv) = 0. (8) 



3 Numerical algorithms to solve relativistic hydrody- 
namics 

In numerical applications, one solves finite difference versions of eqs. (§|-§p. The algorithms 
discussed here are Eulerian, i.e., the quantities -E, M, and R are discretized on a fixed (and 
in our case, Euclidean) grid in the calculational frame (which is conveniently chosen as 
the global rest frame of the system under consideration). They are also explicit, i.e., the 
calculation of quantities at the (discrete) time step n + 1 involves only quantities at previous 
time steps. 

The first step is to treat the 3-divergence operator in eqs. (^|-§) via the method of time- 
step splitting (operator splitting), i.e., for instance eq. (|^) is actually solved by sequentially 
solving dtR+di{Rv^) = 0, i = x,y, z. (No summation over i is implied. It is advantageous to 
use permutations of the sequence {x,y,z) in each new time step.) Therefore, the numerical 
algorithm is only required to solve equations of the type 

dtU + d,{Uv + f) = 0. (9) 
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3.1 The SHASTA algorithm 

We define for further purpose (lower indices denote cell numbers, upper indices the time 
step) 

A, = f/jVi-^;. (10) 

For the SHASTA |Q, one first computes so-called transported and diffused quantities 

Uj = 1 A, - i A,_i + (g^ + g_) f/; + A A/ , (ii) 

where 



— V2=Fei _ n+i/2. 

1 ± (cjii - ej) 

Eq. (|TTD results from a simple graphical picture of mass transport over a fixed grid and 
yields therefore a conservative transport scheme [^]. The source term A/ may depend on 
the values of / in cell j and the neighbouring cells, like in the original explicit SHASTA 
proposed by Boris and Book 



A/ = -QAf;:i'' - /r'^') - Q-ui^"" - fiT) , (13) 



however, we found that computing the source term with the much simpler prescription 

1 
2 



M = -UflT - fl-i'") (14) 

yields results superior to (|13D (see Section 5 and Appendix). The source terms and the 
velocities entering the are computed at half time steps n + 1/2 (half-step method) to 
ensure second order accuracy in time. Although much less time-consuming, simple first 
order full-step methods (which take values at full time steps n only) yield overshoots at 
shock fronts (e.g. as occurring in the one-dimensional shock model studied in [^) and 
should therefore not be used. 

After transporting the quantity U one has to remove the numerical diffusion inherent in 
the transport scheme (^). Let us define the quantity 

A, = - . (15) 

The usual estimate [30| for the diffusion in ([TTl) leads to the antidiffusion fluxes 

^r = ^4, (16) 

which are subtracted from (^) to yield the flnal time-advanced quantities U^j^^ (cf. eq. ([T9| ) 
below). The superscript indicates "explicit", since this antidiffusion was used in the original 
explicit SHASTA However, we found the antidiffusion of the so-called "phoenical" 

SHASTA to yield better results (see Section 5 and Appendix). For this version, the 
antidiffusion flux reads 

< = ^ (A,-i[A,+i-2A,+A,_i]) . (17) 
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Of course, both for the exphcit and the phoenical SHASTA, the antidiffusion step should not 
create new maxima or minima. This is taken into account by replacing the naive antidiffusion 
fluxes ([l^, |l^) by the "flux corrected" antidiffusion fluxes 

= cr ■ max |o, min aAj^i,\Aj\,aAj_i | , (18) 

a = sgn Aj. The final time-advanced quantities read 

rn+l 



m^' = U,-A,+A,.,. (19) 



The variables Uj should fulfill the relativistic constraint E > |M|. A simple way to ensure 
this is discussed at the end of the next subsection. 



3.2 The relativistic HLLE algorithm 



The relativistic HLLE algorithm was originally proposed in p8| and is a Godunov-type 



algorithm |^ . A so-called Godunov algorithm |31[] makes use of the fact that hydrodynamic 
fields are piecewise constant inside each cell and discontinuous at the cell boundaries after 
each time step. The exact time evolution for this particular kind of initial condition is given 
by the decay of all these discontinuities. The decay of a discontinuity is the so-called Riemann 
problem and is in principle analytically solvable (see, for instance, Ref. |^ and references 
therein). A Godunov algorithm simply employs this solution of the Riemann problem for each 
cell boundary to solve the hydrodynamical equations over the whole calculational grid (see, 
for instance, [0]). In turn, Godunov-type algorithms |31|, ^ do not use the exact solution 



of the Riemann problem at a cell boundary but approximate it in a manner consistent with 
the integral conservation laws. 

To be more specific, we consider the cell boundary at x = and the initial distribution 
of the field U at time t = 0, 

with Ui ^ Ur- For times t > the discontinuity at a; = will decay and produce a distribution 

Ui , X < bit 
U{x,t) = { Uir{x) , bit<x<brt (21) 

Ur , X >bj.t . 

bi < and br > are the so-called signal velocities with which the decay proceeds to the 
left or right, respectively. It is obvious that one should consider times t < Ax/|6; ,.|, since 
otherwise the decay of the discontinuity reaches the next cell boundary. If the algorithm 
fulfills the CFL criterion A = At /Ax < 1 this is automatically avoided (since the signal 
velocities are causal). 

In the exact solution to the Riemann problem, Uir{x) is a nontrivial function of x (in- 
volving a simple rarefaction wave, a contact discontinuity and a shock wave, separated by 
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regions of constant flow, see e.g. The relativistic HLLE algorithm approximates this 

function by a constant, which is determined by inserting (|21|) into eq. (^) and integrating 

over a fixed interval [Xmini Xmax]^ ^min < bit, Xmax > bj.t: 

rr - brUr - biUi - FjUr) + FjUi) 

= b;^i ' ^22^ 

where F{U) = Uv + f. Integrating eq. (H) over the (fixed) interval [xmin, 0] or [0, Xmax] and 
using ( p2D yields a value for F{Uir): 

^Ijj . brFm-kFiUr)+bMUr-Ul) 

^ K'^lr) = ^ Z • l^-^j 



-'r 



To determine time-advanced values Uj for each cell in the relativistic HLLE algorithm, 
one first writes @ in the finite difference form 

f/;+^ = U--X - . (24) 

Taking 6*^+1/2 = -^(^7+1/2)5 where the latter is given by eq. (^3]), one obtains a HLLE scheme 
with first order accuracy in time. To achieve second order accuracy (as for the SHASTA), 
half-step updated values for the numerical fluxes G are determined as follows ||28|. First, cell 
interface values at time step n are obtained via replacing the piecewise constant distribution 
of the by a piecewise linear one, 

t/;i^t/;±^^(A,_i,A,), (25) 

where Aj is defined in eq. (0) and 





for 


lA.-il 


< 






A, 


> 


A, 


for 




> 


|A,| 




A. 


> 



S{Aj_i,Aj) = \ Aj for |Aj_i| > \Aj\ , Aj_i • A^- > (26) 
otherwise 



When applying ( |25D one has also to take care of the fact that the interface values must 
not violate the relativistic constraint E > |M|. Then, half-step updated values Uj'^^^'^ are 
estimated by a simple difference form of (|]), 

U;:'^' = Vh-\ {f{U-+) - F{U--)) ■ (27) 

Again, these values should respect E > |M|. Inserting these values on the right-hand side 
of ( p3D we obtain half-step updated values for Gj+1/2, 

— U 
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The transport problem ( P^ is solved if we specify the signal velocities and ^^+1/2 

the cell interface j + 1/2. We note that if the slope function (^) is globally zero, this second 
order scheme is reduced to first order accuracy in time. 

The physical velocities for signals travelling to the right/left is the relativistic addi- 
tion/subtraction of matter velocity and velocity of sound, {v ± Cs)/{1 ± vcs). It was shown 
in 1^ that a convenient choice for the signal velocities is 

{71+1/2 n+1/2 ^ 

^' T ^olr ' 1 , ri+1/2 n+1/2 ( ' 'y^'^) 
l+VCs l + ^;(.+/)_C,^(._(,)„J 

The max-(min-)function accounts for the proper sign of br (bi). Here 

/p"+l/2 ,"+1/2 /pn+1/2 n+1/2 

f = , , (31) 



+ ?7 



V J+ / n+1/2 n+l/2y 

/ ; , X 2 V^(i+1)- J 

;y^+v^) 

nstead of R in the Roe-averages (|3T|, [3^), si: 



In contrast to p8[ we use instead of R in the Roe-averages (|3T| , [3^ ), since for (net) baryon 



free flow i? = 0. For all problems studied r] = 0.5 in equation ( p2D proves to be a safe choice. 
In regions of vacuum, where the denominator in equations ( pT| , |3^ vanishes, we use the 
simple estimate b^j+1/2 = 1 ^j+1/2 ~ which is robust for all cases considered. 

This completes the description of the relativistic HLLE algorithm. Also here, the final 
variables U^^'^ must obey E > |M|. We found that restricting |M| to E for cells violating 
this constraint provides the simplest way to ensure this. For all problems studied here, the 
effect on the global conservation of energy and momentum is negligible. 



3.3 Transformation between CM and local rest frame 

Both algorithms discussed so far require the values of pressure and velocity at various points 
in the transport algorithm. Since p = p{e, n), the set of local rest frame variables e, v, n has 
to be determined from the calculational frame variables E, M, R. According to the definition 
of the latter, this requires in general an inversion of the five equations 

E = (e + ph^-p, (33) 
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M = (e + p)72v, (34) 

R = n-f . (35) 

Nevertheless, it is not necessary to perform a five-dimensional root search First note 

that M and v point into the same direction. Let |M| = M, |v| = v. Now observe that 
M = {E + p)v and also 

e = E-Mv, n = R{l-vY^\ (36) 

Thus, 

M , , 

(37) 



E + p{E-Mv, R{1 - t;2)i/2) • 

For given calculational frame quantities E, M, R, this is a simple fixed-point equation for 
V and can be solved iteratively. Once v is known, e and n can be inferred from (|36D and 
p from the EoS. The different components of v can be obtained from the collinearity of M 
and V. For a simple EoS, one might even find an analytical solution for v However, 
the advantage of the method described here is that it is independent of the EoS (as long 
as p > 0, so that E > M ensures causal matter flow v < 1) and thus (almost) universally 
applicable. 



4 The one— dimensional expansion into vacuum — ana- 
lytical solution 

In this section we discuss the one-dimensional expansion of semi-infinite, (net) baryon-free 
TN and TA matter into the vacuum. The initial condition is 

e(x,0) = I 
v{x,0) = I 

Other hydrodynamic variables follow from the EoS and relations (|33|^5D. The choice v = 1 
in the vacuum is purely conventional, but it guarantees a continuous hydrodynamic solution 
at the boundary to the vacuum, since in the limit of infinite dilution the velocity of matter 
approaches unity^. 



, X < 
, X > 

X < 
X > . 



5 



(38) 
(39) 



4.1 Relation between thermodynamic properties and type of hy- 
drodynamical solution 

Let us clarify the relationship between the thermodynamic properties of matter, i.e., whether 
it is TN or TA, and the type of hydrodynamical solution. Consider one-dimensional hydro- 
dynamic flow. Then one can show |^ that regions of constant flow can only be separated 

^In multi-dimensional applications, due to the isotropy of the vacuum it is not possible to assign it a 
directed, finite velocity. The only possible choice is then v — 0. 
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by a simple wave (continuous solution of the hydrodynamic equations), or by a shock or a 
contact discontinuity (discontinuous solutions of the hydrodynamic equations). This means 
that between the two regions of constant flow in (^, either the initial discontinuity will 
be preserved, or a single simple wave will form to connect these regions, or a more compli- 
cated flow pattern occurs, involving a sequence of simple waves, discontinuities, and constant 
flow. Unless the pressure p(e) vanishes for all e, and matter does not expand at all due to 
the absence of a driving force, the first case can be ruled out. 

Now consider a simple wave, i.e., continuous hydrodynamic flow. This flow respects 
entropy conservation [Q. For isentropic flow, the t— and x— component of the equation of 
motion (|^) can be combined to yield 



d 



f ± c, d 



\dt 1 ± vcs dx ^ 
where the so-called Riemann invariants are 



7^+ = 



(40) 



'^± = y-yo± [ - 
Jen e' 



Cs de^ 



(41) 



Here y = Artanhf is the fluid rapidity. Obviously, the Riemann invariants are constant 

d7^± 



dt 







(42) 



c± 



along the so-called characteristic curves C±{x, t), the positions x±{t) of which are determined 
by integrating 

_ dx± v±Cs . 

w± = — — = — . (43) 

dt l±vcs ^ ^ 

The characteristics are the world-lines of sonic disturbances on the hydrodynamic flow pat- 
tern. For simple waves moving to the right (i.e., f > 0), one can prove that 7^+ = const., 
and it suffices to consider the C_ -characteristic. For simple waves moving to the left (f < 0), 
71- = const., and only the C+-characteristics have to be considered. 

Consider now simple waves moving to the right (left) and a bundle of the corresponding 
non-trivial characteristics (i.e., the C^^-characteristics) in the x — t— plane. For constant time 
t, the change with x of the inverse slope ( ^31) of these characteristics in the x — t— plane is 



dwz, 



dx 



v'{l 



t4(i 



{iTvCsY 

For a simple wave moving to the right (left), we deduce from TZ± = const. 

o\ Co 



dv = t(1 



e + p 



de 



and with 



dc. 



1 d'^p 
2^ 



de 



(44) 



(45) 



(46) 
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and eqs. (E, H3) we obtain 



Since causality requires \wz^\ < 1, we arrive at the important conclusion 

sgn w' = =Fsgn [S e'] . (4J 



For a simple rarefaction wave moving to the right, i.e., where e' < 0, the inverse slope (|43| ) 
of adjacent C_-characteristics has to increase for S > 0, and to decrease for S < 0, and 
stays constant for E = 0. This means that for TN matter the slope of the characteristics will 
decrease, i.e., they "fan out" in the x — t-plane and the rarefaction proceeds undisturbed as a 
simple wave |^3|, i.e., for TN matter the stable hydrodynamic solution is a simple rarefaction 



wave. In turn, the hydrodynamic flow emerging from the initial condition PU] ) will thus 
be a single simple rarefaction wave. 

On the other hand, for TA matter the slope of the characteristics will increase, i.e., they 



are bound to intersect, which is a signal for the formation of shock waves |^3|. Thus, for TA 
matter the stable hydrodynamic solution is a rarefaction shock wave. This also means that 
the initial discontinuity ( p8| , |39D cannot decay. Since an EoS is usually not completely TA, 
but contains regions where matter is also TN (see subsection 4.3), the hydrodynamic flow 
will feature discontinuities in TA regions as well as simple waves in TN regions of the EoS. 

For the special case S = 0, all characteristics have the same slope which indicates the 
existence of a region of constant flow. If E = globally, we would then simply preserve 



the initial situation (|3^, ^) for all times (as discussed above, the same is true for globally 
TA matter). Since the EoS we will finally consider (see subsection 4.3) has TN regions as 
well, only that part of the initial discontinuity will be preserved which corresponds to energy 
densities where E = 0. This means that even if S does not really change sign but only 
vanishes, for the discontinuous initial conditions ( PB| , ^) we will observe a rarefaction shock 
wave in the subsequent hydrodynamical evolution. 

Similar considerations apply for a simple rarefaction wave moving to the left, where the 
C-t--characteristics have to be considered. The result is again that the hydrodynamically 
stable solution is a simple rarefaction wave for TN matter and a rarefaction shock wave for 
TA matter. 

For compressional waves, one infers from ( ^81) that the picture reverses and characteristics 
are bound to intersect for TN matter and "fan out" for TA matter. Thus, the hydrodynam- 
ically stable solution is a compressional shock wave for TN matter and a simple compression 



wave for TA matter. This will be discussed in detail in [R^ 



4.2 Thermodynamically normal matter 

Consider the simple EoS 

p{e) = e , Cg = const. . (49) 

For = 1/3, this is the EoS of an ultrarelativistic ideal gas. One easily checks that 
S = 2c^(l - cl)/{l + cl)e > 0, i.e., this EoS describes TN matter. 
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As shown in the previous subsection, in TN matter the hydrodynamically stable solution 
for the decay of the initial discontinuity ^9]) is a single simple rarefaction wave, where 
= const. [^. For the initial conditions (|38| , |39D we determine the value of 7^+ at 
X = to be IZ^ = and thus 



7^ 

t -- 



For the EoS (|9D this gives 
v{e) = tanht/(e) 



y 



tanh 



Cs de' 
e' + p(e') 



l + c2 



In 



\2cs/{l+cl) 



1 + (e/eo)2'^^/(i+=') ■ 



(50) 



(51) 



This equation nicely shows the defining property of a simple wave, namely that there is a 
unique relationship between the value of the fluid velocity and its thermodynamic state . 
Furthermore, with the help of (^) for C_, we can now calculate for any given e, < e < eo, 
the position x{t; e) at which this value for the energy density occurs for given t, 



x{t; e) 



v{e) - Cs 
l-v(e) Cs 



t . 



(52) 



Eq. (15^ ) has similarity form, i.e., the profile of the rarefaction wave does not change with 
time when plotted as a function of the similarity variable C, = x/t. 

Let us now complete the solution of the hydrodynamic problem. Causality requires that 
the initial conditions ( p8| , pQ] ) remain unchanged for \C^ \ > 1. Thus, we only have to determine 
the solution in the range — 1 < C ^ 1 (the forward light cone). From eqs. (|5l| , ^) we infer 
that the head of the simple rarefaction wave (the point where the rarefaction of matter starts, 
i.e., where the energy density e starts to fall below cq) travels with the velocity — to the 
left. On the other hand, the base of the rarefaction wave (the point where the vacuum ends, 
i.e., where e starts to acquire non- vanishing values) travels with light velocity f = 1 to the 
right. Thus, the energy density as a function of ( can be written as 



e(C) = Cs 



l + Cs 1 + C 



(l+c2)/2c. 



- 1 < C < -c. 

- c, < C < 1 . 



(53) 



The velocity can then be inferred from 
variables from the EoS (EDI) and 



, or simply from (152|) , and the other hydrodynamic 
With the thermodynamic relation dp = s dT = 



(e + p) dT/T and the EoS (^9]) we get the temperature T/Tq = (e/eg 



In Fig. 1 



we show (a) the time evolution of T/Tq, (b) the T/Tq— profile as a function of C (which is 
time-invariant as explained above), and (c) the T™— profile normalized to the initial pressure 
Po for cl = 1/3. 



4.3 Thermodynamically anomalous matter 

Consider the EoS for a (net) baryon-free QGP of gluons, u and d quarks, described by the 



MIT bag model [35|, and that of an ideal gas of massless pions. Match both EoS's via Gibbs 
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equilibrium conditions, i.e., equate pressure and temperature to obtain a first order phase 
transition at temperature Tc with pressure pc = Pq(Tc) = ph(Tc). This yields 

[ r(^)'-B/p, , T>T, 
p{T)=Pc-\ ;M (54) 

Here r = qq/qh is the ratio of degrees of freedom in the QGP and the pion gas {qq = 37 
and (yfji/ = 3 in our case), B is the MIT bag constant, and = [90B/n'^gHir — l)]^^'*. Note 
that B /pc = r ~ 1. As a function of e, 

r (e-4i?)/3 , e>eQ 
P(e) = I Pc , eQ> e> en (55) 

[ e/3 , ej/ > e , 

where eq = 3pc + ^B = (4r — l)pc, = 3pc- Obviously, 

7 e > eg 

, tQ> t> en (56) 
, > e , 

i.e., from eq. (|l]) we infer 

, e > eg 

, eg > e > e/f (57) 
, > e . 

Thus, the EoS (|55|) has TN regions above eg and below e/^-. In the mixed phase S vanishes. 
As discussed in subsection 4.1, for the initial conditions ( p8| , p9| ) this will lead to the formation 
of a rarefaction shock wave, as occurs in the expansion of TA matter. 

If the initial state ( pH]) has an energy density cq < tu-, matter is TN and the expansion 
is exactly the same as discussed in the previous subsection. If, on the other hand, eo > eg, 
the hydrodynamic solution is a simple rarefaction wave in the TN QGP part of the EoS, 
but as soon as the rarefaction reaches the phase transition point e = eg on the wave, one 
enters the mixed phase region of the EoS, and the initial discontinuity will be preserved as 





a rarefaction shock wave p2|, The energy density of matter flowing into this shock 

is eg. Matter emerges from this shock in a hadronic state which will be determined below 
and the rarefaction continues as a simple wave due to that matter being TN. Finally, if 
eg > Co > in the initial state, rarefaction proceeds directly through the rarefaction shock 
from which hadronic matter emerges in a simple rarefaction wave. 

Let us now determine the state of hadronic matter emerging from the rarefaction shock. 
In the rest frame of a (space-like) shock wave, conservation of energy and momentum can 
be written in the form [^, Q 

= T^^ , = T/i , (58) 
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where the subscript i labels quantities in front of the shock. We denote the enthalpy densities 
behind and in front of the shock hj w = e + p, Wi = ei + pi and the corresponding space- 
like components of the hydrodynamic 4-velocity by m = 7f , Ui = 'jiVi. Then, defining a 
generalized volume for (net) baryon-free matter, X = wu'^/wiuf 
equation |45| from eqs. (|58|) [^]: 



HI, one can derive a Taub 



wX-w,-{p-pi){l+X) = . 



Another consequence of the conservation laws (|58D is the relation 



V 



Pi + e 



(59) 



(60) 



for the velocities of matter going into and coming out of the shock wave. 

Suppose the initial state had eo > eg. Then, the rarefaction shock sets in at pi = pc, = 
eg and Wi = 4{pc + B) = Arpc- From the Taub equation (|59|) one derives that the Taub 
adiabat p{X) is given by 



p{X) = 



Pc- 




1-X 



3X -1 



l<X<r 

r <X <4r-l 



(61) 



where the upper line corresponds to final states in the mixed phase and the lower line to 
those in the hadronic phase. Similarly, one can derive a Taub adiabat for the initial state 
being in the mixed phase p4| : 



p{X) = 



1 

3^0 -X 
3X- 1 



1 < X < (3^0 + l)/4 
(3^0 + l)/4 < X < 3^0 



(62) 



where Eq = eo/en- Again, the upper line is for mixed phase final states, the lower for 
hadronic final states. Fig. 2 shows (a) the Taub adiabat ( |6TD and (b) the Taub adiabat ( p^ ) 
for Eq = 5, 10. Also shown is the Chapman- Jouguet (CJ) point on the respective adiabats 
p9| . It is defined as the point where the chord from the initial state (l,Pc) to the final state 
(X, p(X)) has the same slope as the Taub adiabat. This point represents final states where 
the hadronic fluid emerges with the velocity of sound from the shock and where the entropy 
production through the shock is maximized |]39| . 

In principle, the conservation laws (|S3) allow any state on the Taub adiabat as final state 
for the rarefaction shock. This is in contrast to compression shocks in the initial stage of 
a heavy-ion collision where the final compressed state is uniquely determined by the beam 
energy of the colliding nuclei (see e.g. [0). However, mechanical stability of the shock 



||39| , ^ requires that the chord from the initial state (l,Pc) to the final state (X, p(X)) does 
not intersect the Taub adiabat. Therefore, final states below the CJ point are not allowed. 

Furthermore, not all of the remaining final states above the CJ point correspond to the 
long-term hydrodynamic solution, or in other words, lead to a stationary hydrodynamic 
profile. If the final state lies above the CJ point, matter emerges from the shock at a higher 
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density and, consequently, with a smaller velocity than at the CJ point, where it is equal 
to the velocity of sound. However, the subsequent rarefaction proceeds via a simple wave 
where the velocity of matter relative to the wave profile is the sound velocity, see eq. (^31) . 
This means that e.g. for the head of the rarefaction wave (where it is attached to the shock) 
the flow velocity will be larger than that of matter emerging from the shock. Eventually, 
this excess of matter flow leads to a density decrease at the head of the rarefaction wave, 
i.e., in the final state of the shock, and consequently to an increase of the shock strength 
and the velocity of matter emerging from the shock. This proceeds until the latter reaches 
the velocity of sound and matches the velocity of matter on the rarefaction wave. Now the 
solution becomes stationary, i.e., the profile of the rarefaction wave no longer shifts with 
respect to the shock. Thus, in a stationary situation, the final state of the shock will always 
correspond to the CJ point on the Taub adiabat. Remarkably, this process satisfies the 
principle of entropy maximization. 

Let us now complete the solution for the expansion of TA matter into vacuum. We 
consider first the case when the initial energy density of semi-infinite matter eg > eg. 
According to Fig. 3 (a), we can distinguish five regions of flow. Region I extends over 
— oo < X < xa = —Cgt and corresponds to constant flow with e = eo, v = vq = 0. 
Region II starts at xa and is a simple rarefaction wave in the QGP which terminates at 
xb = tiy^eq) — Cs)/(1 — v{eQ)Cs), where (cf. derivation of eq. (|5TD ) 



^ ' 1 + [(e - B)/{eo - 5)]2^=/(i+-?) ■ ^ ^ 

Region III is a region of constant flow, e = eg, v = f (cq) and extends from xb to xc = Vshi, 
Vsh being the velocity of the rarefaction shock wave, which will be determined subsequently. 
Region IV extends from xc to xd = t and represents a simple rarefaction wave in hadronic 
matter. At the head of this wave, e = ecj = Spcj, v = vcj = {vsh + Cs)/(1 + VghCs) (the 



pressure at the CJ point pcj is given below). On this wave, apphcation of (|41|) for 7Z+ = 
yields now 



v{e) = tanhy(e) = tanh 



In I — I + ycj 



(64) 



1 + iecjJ 

where ycj = Artanh vcj is the fluid rapidity of matter emerging from the shock. Finally, 
region V corresponds to vacuum, e = 0, f = 1. 
To determine Vgh we use eq. (|60D, 



(65) 



Cs Pcj + eg 4r - 1 + pcj/p 
where the pressure at the CJ point is given by 



3r - 1 - - (4r - l)/3 

Pcj = Pc / = ■ (66) 

3r-l + 3Jr2-(4r-l)/3 
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In the rest frame of undisturbed matter, the velocity of matter flowing into the shock is 
v{eQ), thus, the shock velocity in this frame is given by 



Vsh 



(67) 



We observe that the solution is again of similarity type. Thus, the energy density as a 
function of the similarity variable ( reads (in the forward light cone) 



e{0=Pc 



1 - Cc 



1-C 



1 (l+c^)/2c. 



1 + C, 



— r + 1 +r — 1 



■ < 



4r 



1 - Cs 1-C 1 + VCJ 
1 + Cs 1 + C l-vcj 



(l+c^)/2cs 



Pc 



-i<C<- 

-Cs<C< 

1 - v{eQ)cs 

Vsh<C<l ■ 



-Cs 

^(cq) 



l-v{eQ)c 

<C <Vsh 



Other hydrodynamic variables can be obtained from the EoS 
The temperature follows from 



(68) 



, eqs. (|63| , |6J), and 



T 



B 



B 



.eg - 
1 



e > eg 

eg > e> en 



(69) 



In Fig. 4 we show (a) the time evolution of T/Tc, (b) the T/Tc— profile as a function of (, 
and (c) the T°°— profile normalized to the critical pressure pc for an initial energy density 
eo = 603|pc- This corresponds to an initial temperature Tq/Tc = 2. Figs. 4 (d-f) shows the 
same for eo = 50 which is only slightly above eg = (4r — l)pc = 48|pc, with an initial 
temperature Tq/Tc — 1.011. 

Let us now construct the solution for the initial state being in the mixed phase. According 
to Fig. 3 (b), we can distinguish three different regions. Region I is constant flow with e = 
eo = EqCh, V = and extends over —oo < x < xa = Vght- Since the velocity of undisturbed 
matter is zero, the (modulus of the) velocity Vsh of the rarefaction shock travelling into that 
region is identical to the velocity Vi of matter streaming into this shock in its rest frame. 
With the help of (|^) this leads to 



where now 



Vsh 



PCJ = Pc 



-Vi 



^PCj/Pc 



3Eo+pcj/Pc 



9E. 



9E^ - lOEo + 1 



9En - 1 



SJ9E^ 



lOEn + 1 



(70) 



(71) 
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and as before Eq = eo/en = eo/3pc- Region II is a simple rarefaction wave in hadronic matter 
and extends from xa to xb = t. Finally, region III corresponds to vacuum, e = 0, v = 1. In 
terms of (, the solution reads 



e(C) = Pc 




1-C l+vcj 
1 + C l-vcj 



(l+c^)/2c. 



Pc 



-l<(<Vsh 
Vsh<C<l ■ 



(72) 



Here, ecj = Spcj with pcj now given by (^, Vsh is given by (|73), and again vcj = 
{vsh + Cs)/(1 + VshCs)- We will not consider numerical solutions for eo in the mixed phase, 
and thus do not show analytical results as in Fig. 4 for the previous case. The reason is that 
the situation is very similar to that in Figs. 4 (d-f) (see also discussion in Section 6) and 
does not exhibit new insights from the analytical or numerical point of view. 

In a sense, the two cases shown in Fig. 4 represent opposite scenarios, both from the 
physical as well as from the numerical point of view. While the first might be the appropriate 
initial condition for a QGP created at RHIC or LHC, the second is more likely to resemble 
initial conditions at the CERN-SPS. The first features a very fast hadronizing rarefaction 
shock, the second a very slow, stationary one. This is illustrated in Fig. 5 where we show 
the shock velocity as given by ( |67| , [70|) as a function of the initial energy density eo- As 



one observes, for a range of initial energy densities around the boundary of mixed phase 
and QGP the shock has small velocities in the rest frame of undisturbed matter. This has 
interesting consequences for the lifetime of the mixed phase in heavy-ion collisions [0, see 
also Section 6. 



5 The one— dimensional expansion into vacuum — nu- 
merical solution 

In this section we compile numerical results for the expansion of semi-infinite matter into 
vacuum. In Fig. 6 we show the expansion of TN matter described by the EoS (^) for the 
relativistic HLLE with A = 0.99 and for the SHASTA with A = 0.4. The grid spacing is taken 
to be Ax = 1 for both cases. Figs. 6 (a,d) show the numerically calculated T/To-profile for 
0, 10, 20,..., 100 time steps, while the corresponding parts (b,c,e,f) show T/Tq and T°°/po as 
a function of ( after 5, 10, 20 and 50 time steps (for the SHASTA also after 100 time steps) 
as compared to the analytical result of Figs. 1 (b,c). The analytical solution is remarkably 
well reproduced, up to effects due to numerical viscosity. The SHASTA shows a tendency 
to produce small fluctuations around the point where the profile is stationary {v ~ Cg). As a 
long-term run has proved, these fluctuations are not instabilities. They can be removed by 
reducing the antidiffusion, see Appendix for details. 

The advantage of plotting quantities as functions of the similarity variable ( becomes 
obvious from Figs. 6 (b,c,e,f): curves for different times should all lie on top of each other. 
Thus, one can monitor the approach to the analytical solution in an elegant and compact 
way. One observes that this approach is rather fast for the HLLE (about 20 time steps) and 



19 



slightly slower for the SHASTA (about 50 time steps). Of course, since the rarefaction wave 
spreads over a growing number of cells (which is linear proportional to the elapsed time), the 
trivial effect of a better resolution on the ^-axis is partly responsible for this approach. In this 
sense, for the same Ax the relativistic HLLE with the larger CFL number A has an advantage 
over the SHASTA: the physical time has advanced further after the same number of time 
steps than in the SHASTA and, consequently, the similarity wave profile extends over a larger 
number of cells. This can also be seen comparing the mutual distance of symbols (time steps 
5, 10, and 20) in Fig. 6 (b,c) with that in (e,f). Taking this effect into account, we conclude 
that the approach of the relativistic HLLE and the SHASTA to the analytical solution takes 
approximately the same physical time t ~ 20 Xhlle^x ~ 50 XshastA'^x ~ 20 Ax. 

Fig. 7 shows the expansion of semi-infinite TA matter for an initial temperature Tq/Tc = 2 
(cf. Figs. 4 (a~c)) for the relativistic HLLE and the SHASTA, respectively. The reproduction 
of the complicated hydrodynamic solution is quite good for the HLLE and satisfactory for 
the SHASTA, although the latter has, after 100 time steps, still problems separating the 
simple rarefaction wave in the hadronic phase from the rarefaction shock and again shows 
wiggles where the profile is stationary. However, as in Fig. 6, the approach to the analytical 
solution takes about the same physical time for both algorithms. 

In Fig. 8 the performance of relativistic HLLE and SHASTA is investigated for the initial 
condition of Figs. 4 (d-f). While the approximation of the analytical solution is quite good 
after 100 time steps for both algorithms, the SHASTA performs globally better for this 
problem than the HLLE. In particular, it approaches the analytical solution after a shorter 
physical time. 

In Fig. 8 (b) it seems that the HLLE algorithm places the rarefaction shock wave at too 
large values of ( in the first few time steps and approaches the correct position only for later 
times. Inspecting Fig. 8 (c), however, reveals that this interpretation is misleading. Here 
the discontinuity is correctly placed, but smeared out due to numerical diffusion. While this 
smearing is symmetric around the true position of the discontinuity in the energy density, it 
appears asymmetric in the temperature, because this quantity is proportional to the fourth 
root of the energy density, cf. eq. (|6^). 

To quantify the deviation from the analytical solution we plot in Fig. 9 the quantity 

d ^ £ dc [r°L(C) - T°1(C)] ' (73) 

as a function of the number of time steps. (The integral is approximated by a Riemann sum, 
/ dx /(x) ~ J2i f{xi){xi+i — Xj)). We also measure the deviation separately in the acausal 
region —2 < C < 1 < C < 3, in order to determine the prediffusion. The corresponding 
part of the deviation (l73| ) is the dotted line in Fig. 9. 

The absolut^ measure of deviation (|73|) confirms the observations made in Figs. 6-8, 
namely that the HLLE is more accurate than the SHASTA, except for the expansion problem 
of Fig. 8. For the problems of Figs. 6, 7 this holds even when comparing the deviation at 
the same physical time. 



^There is no reasonable way to define a relative deviation due to the occurrence of vacuum. 
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Furthermore, for the HLLE the prediffusion is essentially zero up to 50 time steps. At 
this time, the prediffusion has just "leaked" half a cell size (50(1 — X)Ax = Ax/2) into 
the acausal region. The astonishing fact is that, despite the rather disadvantageous A = 

0. 4, the prediffusion in the SHASTA is remarkably small (cf. also the profiles in Figs. 6- 
8). A calculation with the HLLE for the same A = 0.4 gives a much larger deviation 
due to prediffusion than that for the SHASTA. Note that the complexity of the expansion 
solution for TA matter as compared to TN matter has caused to increase d by several 
orders of magnitude in Figs. 9 (b,c) as compared to Fig. 9 (a). However, since d is not a 
relative measure for the deviation from the analytical solution, it can only serve for mutually 
comparing the performance of different algorithms. 

The above investigation has shown that the algorithms are well able to handle the vacuum 
and rarefaction shock waves in the expansion of TA matter. However, they need a certain 
number of time steps before they approach the analytical solution. If one wants to limit the 
initial deviation to a given physical time interval one has to choose a sufficiently small time 
step width At, or, for fixed A, a sufficiently fine spatial resolution Ax. For instance, if the 
numerical reproduction of the hydrodynamical solution is taken to be sufficiently accurate 
after 50 time steps and one wants this accuracy to be achieved after a physical time of 1 
fm, one is bound to choose At = lfm/50 = 0.02 fm, or Ax = 0.02 fm/A. A system with 
longitudinal extension L must then be resolved over 50AL/fm grid points. This requires 
rather large grids in order to study dynamical problems over longer periods of time. 

Let us finally mention that in the HLLE runs we do not use the physical velocity of 
sound in the mixed phase (c^ = 0) in the signal velocity estimates (|2^, ^). We found that 
doing so prevents matter transport across the rarefaction shock if the velocity of the latter 
is small and leads to creation of vacuum between rarefaction shock and hadronic simple 
wave. Rather, we globally take = 1/3 in the signal velocity estimates which is safe but 
also increases the numerical dissipation. Furthermore, in the actual SHASTA propagation, 
we take v = for vacuum cells, as is also necessary for multidimensional applications, cf. 
footnote 6. We checked that the difference to the choice w = 1 as in (^) is negligible anyway. 

6 Expansion in finite systems 

In this section we discuss the one-dimensional expansion into the vacuum for finite systems, 

1. e., we start with the initial condition 

{—1 , — oo < X < —R 
, -R<x<R (75) 
1 , i? < X < oo , 

where 2R is the size of the system {R is the "radius"). Let us briefly sketch the subsequent 
hydrodynamical evolution. In Section 4 we have seen that a simple rarefaction wave travels 
into undisturbed TN matter with velocity Vrare = Cs- For an initial state in the mixed phase. 
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the rarefaction is accomplished by a shock with velocity Vrare = Vsh as given by (|7D|). Thus, 
for times t < R/vrare, i-e., as long as the rarefaction has not reached the center of the system 
at X = 0, the evolution of the system for x > is exactly the same as described in Section 
4, and for x < simply the reflected version of that. In this stage, the solution is still of 
similarity type. However, as soon as the two rarefaction waves reach the center and start to 
overlap, the solution becomes more complicated and is no longer of similarity type. 

Nevertheless, as was shown by Landau IQ, if matter has an EoS of the form (^), an 
analytical solution still exists. Its derivation is lengthy and shall not be repeated here. 
The general structure, however, is the following: in the center of the system where the (in 
this case, simple) rarefaction waves overlap, the hydrodynamic solution is determined by 
a potential xiT,y). As matter expands, the overlap region grows, but since it can do so 
only causally, outside this region the hydrodynamic solution remains a simple rarefaction 
wave. The boundary between the two regions is determined by smoothly matching the 
hydrodynamic flow. 

A practical way to generate a semi-analytic solution to this problem in the whole forward 
light cone is via the method of characteristics |21|]. The basic idea is to solve the system of 



four characteristic equations (^, ^o|) which are ordinary differential equations. This method 
was described in detail in ||21[] and shall not be repeated here. We will use it to compare 
with the results produced by the algorithms of Section 3. This method, however, is also 



restricted to TN matter, since the characteristics will intersect at shock discontinuities [|3 
as occur in the expansion of TA matter. Then, the initial value problem posed by (|42| , |43| ) 
is no longer unique and thus not solvable. 

In Figs. 10, 11 we show the results for the expansion of TN matter with = 1/3 for 
the relativistic HLLE and the SHASTA in comparison to the semi-analytic result generated 
by the method of characteristics |^1|. We again choose A = 0.99 for the HLLE and A = 



0.4 for the SHASTA. Of special interest is the dependence of the results on the spatial 
resolution: Figs. 10, 11 (a,b) show the temperature and (c,d) the energy density for Ax = 
0.1 i? and Ax = 0.025 R, respectively. One observes that in the first case both quantities are 
overestimated in the central region at intermediate times and also that viscosity effects are 
more prominent. 

This is explained as follows: as the results of the previous section have indicated, it 
takes about 20 time steps for the HLLE and of the order of 50 for the SHASTA to correctly 
approximate the analytical shape of a simple rarefaction wave. This criterion is not fulfilled 
in Figs. 10, 11 (a,c), since the simple rarefaction wave reaches the center after about 14 time 
steps in Figs. 10 (a,c) and about 30 time steps in Figs. 11 (a,c). The two rarefaction waves 
then start to overlap before they were able to assume their correct shape. Subsequently, 
the numerics fail to reproduce the correct behaviour in the overlap region, leading to the 
overestimate of temperature and energy density in that region. 

On the other hand, in Figs. 10 (b,d) the simple rarefaction wave travelled for about 
60 time steps (HLLE) and in Figs. 11 (b,d) for about 130 time steps (SHASTA) before 
overlapping with its partner. It thus had sufficient time to establish the correct shape. 
Consequently, the analytical values for quantities in the overlap region are much better 
approximated by the numerics. As in the expansion of semi-infinite matter, the relativistic 
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HLLE performs better than the SHASTA. We note that in all cases energy and momentum 
is globally conserved on the level of the calculational accuracy (i.e., 10~^^). 

From the physical point of view we would like to draw attention to the phenomenon that 
for fixed time t in the global rest frame of the system, matter at small |x| in the overlap 
region of the two rarefaction waves is colder and less dense than for larger This is a 
purely relativistic effect: since the matter velocity is a monotonously increasing function of 
and since time in the local rest frame of fast matter is dilated as compared to that of 
slow matter, at a fixed global t the rarefaction for cells at small |x| has already proceeded 
further than for cells at larger \x\. This leads to a localized wave of denser material travelling 
outward which is best seen in the energy density profiles. Figs. 10, 11 (c,d). 

In Fig. 12 we show the expansion of matter described by the EoS ( p5D for the relativistic 
HLLE and in Fig. 13 for the SHASTA. Part (a) and (c) show temperature and energy density, 
respectively, for cq = 603|pc, while part (b) and (d) show these quantities for eo = 50 Pc- We 
choose Ax = 0.025 i?, which, according to the results of Section 5, allows for a reasonable 
reproduction of the analytic similarity solution before overlap occurs. In order to better 
distinguish the profiles for different times we show them alternatingly as full and dotted 
lines. As long as the rarefaction waves do not overlap, the codes reproduce the features of 
the expansion of semi-infinite TA matter discussed in Sections 4, 5. The SHASTA again tends 
to produce small ripples at the stationary point of the QGP rarefaction profile for the case 
To/Tc = 2, see Figs. 13 (a,c), which can again be suppressed by reducing the antidiffusion. 
This also removes the instabilities on the hadronic rarefaction wave visible in Fig. 13 (a). 
Global energy-momentum conservation is again fulfilled within calculational accuracy. 

Figs. 12, 13 give a clear picture of the difference in the physical scenario if we initialize a 
very hot system {Tq/T^ ~ 2) or a system near the phase transition temperature {Tq/Tc ~ 1). 
Suppose eo > cq. Then, the simple rarefaction wave accelerates QGP matter until the 
rarefaction shock sets in when e = eg. The final hadronic state is always given by the CJ 
point (|66|) . The important point is that initial and final state of the shock are independent 
of eo. Thus, the velocity of the shock in the rest frame of matter moving into the shock is 
independent of eo. 

However, for large eo, the acceleration of matter in front of the shock is large, so that 
the latter is explosively driven outwards by the expanding matter. On the other hand, for 
eo not much larger than eg, the acceleration of QGP in front of the shock is small, and 
consequently, the velocity of the latter is small or even negative (cf. Fig. 5). The system 
"burns out" from the inside and seems to "evaporate" hadrons. In both cases, the energy 
density in the interior decreases until e drops below en and hadronization proceeds also from 
the inside. As explained above, for a fixed time in the global rest frame of the system, this 
happens in general earlier for cells at small 

Let us mention that for eo in the mixed phase, there is no pre-acceleration by a simple 
rarefaction wave and the rarefaction shock "burns up" mixed matter with a velocity Vsh 
(given by eq. ([70D) that is the larger the smaller eo is. For eo = ej/, this velocity reaches —Cg 
(cf. Fig. 5) and the rarefaction shock degenerates to a simple rarefaction wave. There is no 
hadronization from the inside in this case since mixed phase matter is at rest and simply 
vanishes as soon as the two rarefaction shocks coming in from right and left meet at the 
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center. 

This poses the question whether a remnant of the QGP survives longer if the system 
is initially very hot or rather, if it starts its evolution near the phase transition. In more 
physical terms one could ask what is the optimum beam energy for heavy-ion collisions to 
prepare a quark-gluon system with a maximum lifetime |47|. In 1 47] it was shown in a three- 
dimensional calculation (employing cylindrical symmetry in two space dimensions) that the 
mixed phase survives longer if the system is initialized in this phase than in the hot QGP 
phase. This behaviour was explained by the existence of a so-called "softest point", i.e., a 
minimum m p/e as a function of e, in the mixed phase region of the EoS. The reduction of 
the pressure relative to the energy density was held responsible for the observed small flow 
velocities of mixed phase matter. Thus, in the rest frame of the system, mixed phase matter 
does not expand and is only slowly consumed by "evaporating" hadronic matter with high 
velocity from its surface (see Fig. 2 of Ref. p7| ). 

From our above investigations we can sustain this picture with a more detailed under- 
standing of the slow burning of mixed phase matter and the "evaporation" of fast hadronic 
matter as being due to a rarefaction shock wave. To answer the question about the maximum 
lifetime of the mixed phase we show in Figs. 14 (a,b) isotherms for a temperature Tc — 5 (or 
equivalently, curves of constant e = en — 5', 5,5' arbitrarily small), as calculated with the 
relativistic HLLE (each dot corresponds to a single cell with the respective temperature). 
One observes that for fixed time, the "burning" of the mixed phase proceeds exclusively 
from the outside for eo near eg while for larger eo hadronization can also occur from the 
inside. This has the consequence that, in agreement with the conclusions of [J^, at least at 
X = (around midrapidity) the lifetime of the mixed phase is shortened by the rarefaction 
following an "explosive" expansion scenario and is prolonged in a "slow burning" scenario. 
(The apparent long lifetime of the mixed phase at finite x as observed for eo > 100 is 
again the mentioned relativistic time dilation effect: for fixed t the proper time of the cor- 
responding fluid elements has simply not very much advanced due to their high velocity in 
the calculational frame; therefore, they are still hotter.) 

However, this is only true up to a certain initial energy density: if eo is very large, it 
also takes longer to cool the interior and the lifetime starts to grow again. To be specific 
we plot in Fig. 14 (c) the time tfmai where the mixed phase ceases to exist at a; = 0. For 

^ ^0 ^ ^H, mixed phase matter is consumed by a rarefaction shock with constant velocity 
Vsh given by ([701). Thus, the (Tc — 5)-isotherm can be analytically calculated: it is simply a 
triangle with its tip located at x = and t = t final I R = '^l\vsh\- Since \vsh\ decreases when 
eo approaches eg (cf. Fig. 5), the lifetime increases until eo = eg. For larger eo there is no 
longer an analytical solution and we have to read off the lifetime from Fig. 14 (a) (diamonds 
in Fig. 14 (c)). However, for very large eo we expect the hydrodynamical behaviour at x = 
to resemble that of Bjorken's solution [^, i.e., t final ~ ■ Such a behaviour is given by 



the dotted line in Fig. 14 (c) and indeed confirms our expectation. 

As one observes, there is a local maximum in the lifetime of the mixed phase, as predicted 
|7|, but asymptotically the lifetime grows indefinitely proportional to e^J'^. In the realistic 



m 



three-dimensional case, the additional spatial dimensions will probably only affect the global 
time scale of cooling, and thus we expect this picture to remain qualitatively correct. 
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7 Conclusions 



In this work we have presented two numerical schemes, the SHASTA and the relativistic 
HLLE, to solve the equations of ideal relativistic hydrodynamics. These algorithms can be 
applied to the simulation of heavy-ion collisions in realistic 3+1-dimensional space-time 
geometry, i.e., without assuming specific symmetries to simplify the solution. We anticipate 
that they will find application to analyze recent data in the BEVALAC energy range 0, 
in order to understand to what extent ideal relativistic hydrodynamics describes heavy-ion 
collisions from the initial compression stage to the final freeze-out. Furthermore, they are 
useful as underlying transport algorithms for multi-fluid-dynamics W^ which is applicable 



in the AGS to SPS energy range |]I4[. Finally, they are applied to study the dynamics of 
spatial inhomogeneities created in the central region of ultrarelativistic heavy-ion collisions 
at RHIC energies 

We have in detail presented the analytical solution for the one-dimensional expansion 
into vacuum, both for thermodynamically normal as well as thermodynamically anomalous 
matter. We have shown that the algorithms reproduce the analytical solution rather well 
after a certain number of time steps. Choosing a sufficiently small spatial grid size Ax, 
the physical time until the analytical solution is reproduced can in principle be made ar- 
bitrarily small. In realistic situations, this will be ultimately limited only by the available 
computational power. 

For the one-dimensional expansion of finite systems we have investigated the lifetime of 
a mixed phase of pions and QGP. As expected from scaling hydrodynamics, for sufficiently 
large initial energy density eo this lifetime grows ~ Bq^^. However, the lifetime has also a 
local maximum near eg, i.e., around the phase boundary between QGP and mixed phase 
matter which opens the opportunity to study a long-lived system of deconfined matter not 
only in the ultrarelativistic energy domain but probably also by tuning the beam energy of 



presently available accelerators Let us nevertheless note that in order to detect this 

long-lived system via e.g. its electromagnetic radiation, not its lifetime at x = 0, nor the size 
of its space-time volume is most decisive, but rather its initial temperature: the number of 
emitted photons or dileptons is only linear proportional to the space-time volume, but grows 
proportional to the fourth power of the temperature. A hot system is more efficiently "out- 
shining" other sources of electromagnetic radiation than a large, but relatively cool system. 
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A Appendix 



This appendix serves to illustrate the effect of modifications of the SHASTA. We first demon- 
strate why the phoenical SHASTA (defined by the antidiffusion fluxes (|T7D) with the simpli- 
fied source treatment is preferable to the phoenical SHASTA with the original sources 
(p!3|) and the explicit SHASTA (with the antidiffusion fluxes (p!6D ) with either original or 
simplified sources. In Fig. 15 we show results for the expansion of (a-c) TN matter and (d- 
f) TA matter with eo = 50 Pc as calculated with these alternative versions. As one observes 
in Fig. 15 (a), the phoenical SHASTA produces a "terrace" near the stationary point of the 
simple wave profile. The explicit SHASTA versions perform satisfactorily for this expansion 
problem (Figs. 15 (b,c)). However, all three versions fail to produce the correct profiles for 
the expansion of TA matter, see Figs. 15 (d-f). 

We furthermore show how one can suppress "wiggles" on the stationary point of the rare- 
faction wave profile as produced by the SHASTA version used in the main part of this work, 
cf. Figs. 16 (a,b). To this end, we simply reduce the antidiffusion fluxes {\i7\) by replacing 
the first factor 1/8 by a smaller number, e.g. 1/10. The result is shown in Figs. 16 (c,d). 
The effective reduction of the antidiffusion by 20% causes the wiggles to disappear, but, of 
course, it also increases the numerical diffusion, especially visible in the larger prediffusion 
into the vacuum. Nevertheless, such a reduction of the antidiffusion fluxes proves to be 



crucial to stabilize the transport algorithm in two space dimensions for large times [27 
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Figure Captions: 



Fig. 1: Analytic solution for the expansion of TN matter with EoS (^). (a) Tempera- 
ture normalized to the initial temperature as function of the spatial variable x for times 
t = 0, 10, 20, 100. Since the problem is of the similarity type, x and t carry no units 
and can be scaled with an arbitrary, constant factor, (b) T/Tq, and (c) energy density 
normalized to the initial pressure po as functions of the similarity variable (. 

Fig. 2: (a) The Taub adiabat (|6lD . (b) The Taub adiabat (0) for £"0 = 5 (solid line) 
and Eq = 10 (dotted hne). The corresponding Chapman- Jouguet points and the chord 
connecting them with the center {X,p) = (l,Pc) of the adiabat are also shown. 

Fig. 3: Schematic rarefaction profile in TA matter for initial states (a) eo > eg and (b) 
eg > Co > See text for a detailed discussion of the different regions. 

Fig. 4: Analytic solution for the expansion of TA matter with EoS (^) for (a-c) eg = 603|pc 
and (d-f) eo = 50 pc- (a,d) T/Tc as function of x for times t = 0, 10,20, 100, (b,e) T/Tc, 
and (c,f) T^^/pc as functions of the similarity variable (. 

Fig. 5: Velocity of the rarefaction shock in the rest frame of undisturbed matter as a 
function of the initial energy density eo (in units of pc). 

Fig. 6: Numerical solution for the expansion of TN matter with EoS (^), generated with 
(a-c) the relativistic HLLE (Ax = 1, A = 0.99) and (d-f) the SHASTA (Ax = 1, A = 0.4). 
(a,d) T/Tq as function of x for 0, 10, 20,..., 100 time steps, (b,e) T/Tq, and (d,f) T°°/po 
as functions of C after 5 (stars), 10 (squares), 20 (diamonds), 50 time steps (dotted line) in 
comparison to the analytical solution (solid line). The dashed line in (e,f) corresponds to 
the profile after 100 time steps. 

Fig. 7: Numerical solution for the expansion of TA matter with EoS (|55|) for eo = 603|pc; 
generated with (a-c) the relativistic HLLE (Ax = 1, A = 0.99) and (d-f) the SHASTA 
(Ax = 1, A = 0.4). (a,d) T/Tc as function of x for 0, 10, 20,..., 100 time steps, (b,e) T/T^ 
and (c,f) T^^ /pc as functions of C after 5 (stars), 10 (squares), 20 (diamonds), 50 (dotted 
line), 100 time steps (dashed line) in comparison to the analytical solution (solid line). 

Fig. 8: As in Fig. 7 for eo = 50 pc- 

Fig. 9: (a-c) Deviation ([73D for the expansion problems of Figs. 6-8. Solid lines are the 
total deviation, dotted lines the deviation outside the causal light cone; circles correspond 
to the relativistic HLLE, squares to the SHASTA. 

Fig. 10: Expansion of a finite system of TN matter with EoS (P5|), calculated with the 
relativistic HLLE (A = 0.99). (a,b) T/Tq- and (c,d) T°°/]9o^profiles for times t = nXR, 
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n = 0, 1,2, 10. For (a,c) Ax = 0.1 i?, for (b,d) Ax = 0.025 R. In (a,c) each dot corre- 
sponds to a cell (in (b,d) a similar presentation is not possible since individual dots would 
not be resolved on the scale of the picture; we therefore choose to draw the profiles as dotted 
lines) . 

Fig. 11: As in Fig. 10, calculated with the SHASTA (A = 0.4). 

Fig. 12: Expansion of a finite system of TA matter with EoS (|5BD, calculated with the 
relativistic HLLE (Ax = 0.025 i?, A = 0.99). (a,b) T/T^- and (c,d) T^Vp^-profiles for 
times t = n\R, n = 0,1, 2, ...,10 (alternatingly shown as full and dotted lines). For (a,c) 
eo = 603|pc, for (b,d) eo = 50 pc- 

Fig. 13: As in Fig. 12, calculated with the SHASTA (A = 0.4). 

Fig. 14: (a) (T, - (5)-isotherms in the x - t-plane for eo = 50, 52, 60, 100, 200, 603|pc (from 
inward to outward lying curves), calculated with the relativistic HLLE (Ax = 0.025 i?, A = 
0.99). Each dot corresponds to a cell with the respective temperature, (b) As in (a) on a 
larger scale, in order to show the complete shape of the isotherm for eo = 603|pc- (c) The 
lifetime of the mixed phase (i.e. the intersection of the isotherms in Fig. 14 (a) with the 
t-axis) as a function of eo. The sohd curve (corresponding to eg > eo > e^) is calculated 
analytically. The diamonds are the values as read off in Fig. 14 (a). The dotted line indi- 
cates the eo^^-behaviour expected for asymptotically large eo from Bjorken's scaling solution. 

Fig. 15: Temperature profiles for the expansion of (a-c) TN matter and (d-f) TA mat- 
ter with eo = 50pc for 0, 10, 20,..., 100 time steps; (a,d) phoenical SHASTA with original 
source treatment, (b,e) explicit SHASTA with original source treatment and (c,f) with sim- 
plified source treatment. For all cases Ax = 1, A = 0.4. 

Fig. 16: Temperature profiles after and 100 time steps for the expansion of (a,c) TN 
matter and (b,d) TA matter with eo = 50 Pc, calculated with the phoenical SHASTA with 
(a,b) regular antidiffusion and (c,d) antidiffusion reduced by 20% (Ax = 1, A = 0.4). 
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